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Recently Han and Heary proposed an approach to steady-state quantum transport through meso- 
scopic structures, which maps the non-equilibrium problem onto a family of auxiliary quantum 
impurity systems subject to imaginary voltages. We employ continuous-time quantum Monte-Carlo 
solvers to calculate accurate imaginary time data for the auxiliary models. The spectral function is 
obtained from a maximum entropy analytical continuation in both Matsubara frequency and com- 
plexified voltage. To enable the analytical continuation we construct a kernel which is compatible 
with the analytical structure of the theory. While it remains a formidable task to extract reliable 
spectral functions from this unbiased procedure, particularly for large voltages, our results indicate 
that the method in principle yields results in agreement with those obtained by other methods. 

PACS numbers: 72.10.Bg, 73.63.Kv 



I. INTRODUCTION 



The calculation of steady-state transport properties 
of open quantum systems such as quantum dots is a 
challenging and unsolved problem. Pcrturbative meth- 
ods [TH3] may be used to study the weak correlation 
regime, but they fail to provide a reliable description of 
the competition between Kondo- and Coulomb-blockade 
physics in strongly interacting dots [4]. To avoid these 
limitations of conventional perturbation theory, various 
non-perturbative numerical approaches have been devel- 
oped. Time-dependent density-matrix renormalization 
group (tDMRG) calculations [H [6] and real-time Monte 
Carlo (RT-MC) approaches [7HTU] try to compute the 
relaxation into the interacting steady state after some 
switching of parameters, such as voltage bias or interac- 
tion. While the short-time transients can be very accu- 
rately captured with these methods |11) , the approach to 
the steady-state may occur on rather long, in the worst 
case exponentially large times scales. Due to finite-size 
effects in the tDMRG and an exponentially growing sign 
problem with increasing time in RT-MC, the access to 
long times is severely limited in both approaches. Fur- 
thermore, the tDMRG is performed for a finite, closed 
system; whether a relaxation to a reasonable approxima- 
tion of the interacting steady-state is guaranteed for some 
intermediate time scale much smaller than Poincare's re- 
currence time is not obvious. This latter problems may 
be avoided by numerical renormalization group (NRG) 
[P2] and functional renormalization group (fRG) calcula- 
tions [P3HT8] . which attempt a direct description of the 
non-equilibrium steady state. However, the former in- 
troduces an artificial discretization and truncation of the 
spectrum of the Hamiltonian, which can lead to artifacts 
in the time evolution. The fRG, on the other hand, is 
again perturbative in nature, and experience up to now 
shows that it works best in the extreme non-equilibrium 



limit p7] . 

None of the methods developed so far is able to pro- 
vide a complete and reliable description of the model in 
all parameter regimes. More importantly, the most inter- 
esting regime, where all relevant energy scales - voltage, 
temperature, magnetic field etc. - are of the same order 
as the relevant low-energy scale of the model, is usually 
the one which is not accessible. Therefore, the devel- 
opment of new or improved simulation approaches is a 
worthwhile and important task. 

Recently, a new and rather unconventional approach to 
calculate the steady-state transport through interacting 
quantum dots or similar structures was proposed by Han 
and Heary 19J . Their formalism, which is based on Her- 
shficld's density operator |20| . maps the non-equilibrium 
steady-state of the interacting model onto an infinite set 
of auxiliary equilibrium systems, each characterized by 
some complex voltage. The appealing feature of this ap- 
proach is that powerful methods exist for the numerical 
solution of equilibrium models. The complexification of 
the voltage bias, however, introduces a formidable new 
problem in the form of an analytical continuation in the 
voltage on top of the already challenging analytical con- 
tinuation from Matsubara frequencies to real frequencies. 
In Ref. |19j this double analytical continuation was per- 
formed using a phenomenological formula based on gen- 
eral structures of the self-energy found in second order 
perturbation theory. 

The purpose of this study is to explore to what extent 
an unbiased numerical implementation of the method by 
Han and Heary is feasible. We will address two issues: (i) 
the use of recently developed, accurate continuous-time 
quantum Monte-Carlo (CT-QMC) algorithms to simu- 
late quantum impurity models as solvers for the effec- 
tive equilibrium impurity problems with complex volt- 
age bias; and (ii) the analytical continuation of Matsub- 
ara frequency data via some Maximum Entropy method. 
In particular, we will compare the performance of the 
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weak-coupling [2T] and hybridization expansion [55] al- 
gorithms and propose a kernel for the Maximum Entropy 
(ME) procedure which is compatible with the analytical 
properties of the Green function. 

The paper is organized as follows. Section |TT| describes 
the imaginary-time approach to steady state transport 
by Han and Heary. A brief introduction to the CT-QMC 
for equilibrium problems and their suitability for mod- 
els with complex voltage bias follows in section [TTT| Sec- 
tion |IV B| is devoted to the issue of analytical continua- 
tion in the voltage and frequency domain and presents 
some results for equilibrium and non-equilibrium situa- 
tions. We will finish the paper with a conclusion and 
outlook in section [Vl 



II. IMAGINARY-TIME FORMULATION OF 
STEADY-STATE TRANSPORT 

We briefly review the imaginary-time formulation of 
steady-state transport through an interacting quantum 
dot proposed by Han and Heary [19] . which is based on 
the work of Hcrshfield 1201. 



A. Physical Model 

We consider a spin-degenerate, single-level quantum 
dot attached to two non-interacting fcrmionic leads. This 
system can be described by the Single-Impurity Anderson 
Model with Hamiltonian (e = H = 1) 
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where a = — 1 and a = +1 label the left and right reser- 
voirs, respectively. The index fc denotes the wave-vector 
of the lead states and a the spin quantum number. A 
gate voltage Vq may be applied to shift the dot energy 
level position relative to the particle-hole symmetric con- 
figuration Vg = 0. 

To keep things simple, we assume a fc-independent hy- 
bridization V a ka = V/y/2 and consider the wide-band 
limit for the dispersion of the leads. We then end up with 
a bare level broadening r = T L + T R , T a = ir\V\ 2 N F /2, 
where Np denotes the density of states of the leads at 
the Fermi energy. 

In the case of non-equilibrium steady-state transport, 
the leads are supposed to be unaffected by the cur- 
rent flowing through the dot and characterized by free 
Fermion correlators 



with fp(x) = (e l3x + l)^ 1 the Fermi distribution function 
for inverse temperature /3 and [L a the value of the chemi- 
cal potential for lead a. We restrict ourselves to the case 
where the inverse temperatures of the left and right lead 
are the same, /?/, = /3r = /?, and symmetrically applied 
voltage bias, hl = — M-R- The bias voltage is denoted by 
$ = fi L - n R . 



B. The Y-Operator 

In Ref. [501 Hershfield introduced a Hcrmitian operator 
Y by means of which the non-equilibrium, steady-state 
expectation value of a local observable A may be written 

as 



Tve -/3(H~<S>Y) ■ 



(5) 



The above expectation value is of the form (A) = 
Tr pA/Tr p, and hence resembles the equilibrium expres- 
sion. Under certain assumptions involving the non-trivial 
exchange of limiting procedures, the operator Y can be 
expressed as 



Y = y ^ 

/ j 2 ether' ako 
aka 



(6) 



where the scattering states ip a k<r are related to the 
bare conduction states c a ka by the second-quantized 
Lippmann-Schwinger equation |23j 
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The Liouvillians are defined as C = [H, ■] and Ly — 
[H v ,-], with H v = Y, a ka( v akacl ka d a +h.c.) the hy- 
bridization part of the Hamiltonian. The "•" denotes the 
operators after C, and the fraction in Eq. Q denotes 
the corresponding geometric series in C, i.e. a series of 
iterated commutators with H. 

For U 7^ it is impossible to calculate an explicit ex- 
pression for the Y-operator. More importantly, although 
H — $F looks like an effective Hamiltonian for the sys- 
tem, it cannot be used to define a consistent description 
of imaginary-time and real-time dynamics. The real-time 
dynamics is always controlled by H alone, but H and 
H — <&y will in general have a different spectrum. There- 
fore, the analytically continued imaginary-time dynamics 
does not reproduce the real-time dynamics. 



C. Imaginary Voltages 

Since H — $Y does not yield the correct real-time dy- 
namics, Han and Heary 19J introduce an additional trick. 
Starting with a fully established non-interacting steady- 
state ensemble at time t = 0, the fully interacting steady 
state is formally reached by propagating the system to 
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t = +00. In a path integral representation the expecta- 
tion value for an observable A becomes 





(8) 

Here, the average (-)o is performed using Eq. ^ with 
H —> Hq and Y Yd, where Yq can be explicitly con- 
structed using non-interacting scattering states. It was 
argued in Ref. [THlthat the time evolution via H maps the 
non-interacting scattering states to the interacting ones 
and the Lagrangian for the real-time evolution reads 



(9) 



correct time evolution. 

To obtain a Matsubara-like theory, the fields ip are 
now Wick rotated, tp(t) — > ip(—iT). However, under the 
replacement t — > —it the exponential factor becomes 
e Q * T / 2 , which means that it diverges as r — > 00 and de- 
cays as t — y — 00. To circumvent this problem, Han and 
Heary introduce a second analytic continuation to ensure 
Matsubara's periodic boundary conditions and thereby 
obtain a well-defined effective equilibrium system. This 
is achieved by complexifying the voltage occurring in the 
extra time evolution rate according to $ — > itp m , ra£Z. 
For the particular choice cp m — 4irm/f3 the Matsubara 
boundary conditions are conserved |19j . 



Aiming at a description which yields e lH ^ _t ) as real- 
time evolution operator for t — » t' and c~( r as 
imaginary-time evolution operator for — ir — > —'it', the 
Lagrangian is reexpressed with respect to the spectrum 
of H — $F, e a ka = £ a k<j — a$/2. Statistical expectation 
values take a form analogous to equilibrium expectation 
values, with a uniform Fermi level e a k a — 0. Due to the 
discrepancy between H and H — <&Y , the real-time La- 
grangian transforms to L(t) = Y, a k<? ^Ikai^i^t-Saka- 
a&/2)ip aka (t), so the effective Fermi levels of left and 
right leads have different time evolution rates. These 
rates can be factored out as time-dependent phase factors 
of the Grassmann fields by introducing new field variables 
^Pakcr(t) = e lam ^ 2 tp a ka(t). The extra time evolution rate 
is generated by idt acting on the phase factor, and thus 
L (t) = T,aka^lka( t )( [d t ~ describes the 



D. Effective Action 

The final result of these manipulations is that both 
the Lagrangian and the fields now have their time evo- 
lution with respect to the effective equilibrium Hamil- 
tonian K = H — ($ — i(p m )Y. In a perturbative ex- 
pansion around the non-interacting limit, one may then 
switch to the interaction picture with respect to the non- 
interacting effective Hamiltonian Kq = H^ — (<fr — 'vp m )YQ. 
As before, Yb is Hershfield's boundary condition operator 
for the corresponding fully established non-interacting 
steady state, for which an explicit expression can be 
given. 

We may now proceed along the usual lines and inte- 
grate out the conduction electron degrees of freedom to 
obtain an effective action 



S, 



eff : 



E 



/ f dr dT'dUr')G^(T',T)d a (T) + U ( 
J Jo Jo 



dr 4(rK(r) 



d\(r)d f (T) 



(10) 



for the electrons on the dot. As we are by construction 
in the stationary state, the bare dot Green's function 
Goa( T ' j T ) appearing in the quadratic term in the action 
( 10 1 depends on the time difference only. We therefore 



may perform a Fourier transform to fermionic Matsubara 
frequencies and find the form |19) 



Go,r 



E 



1/2 



a=±l LUJ n 



f(i^ m -$)-e d + ir 



(a) 



, (11) 



with G ,„ m := G (i<y2 m ,kj„), T^l := rsgn(a;„ - a(p m /2), 
and e d = V G . 

The desired Green's function for the stationary state 
of the interacting system is finally obtained by solving 
the quantum impurity problem for each iip m , m € Z, 
performing the analytical continuation icp m z v and 
evaluating the resulting expression at the physical voltage 



Although the preceding discussion seems to be based 
on simple manipulations of the functional integral, one 
has to show formally the equivalence of the complexified 
auxiliary equilibrium time-evolution based on the action 



( 10 ) and the actual physical time evolution with respect 
to H as given by ^ after the analytical continuation 
iip m — > $ in the former. Up to now such a formal proof 
is still lacking, only an argument based on the inspection 
of the contributions to perturbation expansion has been 
put forward [19]. It is therefore interesting to see if an 
unbiased numerical implementation of this formalism is 
possible and produces physically meaningful results. 
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III. CONTINUOUS-TIME QUANTUM MONTE 
CARLO 



In order to compute the self-energy from action 



(10 1 as a function of Matsubara frequency we employ 
continuous-time Monte Carlo (CT-QMC) solvers. The 
continuous-time Monte Carlo technique in the weak- 
coupling |21j and hybridization expansion |22| formula- 
tion has been discussed in considerable detail in the lit- 
erature and we will present here merely a short summary 



of the formalism. The idea is to expand the partition 
function Z = Tr [e~@ H ] into a series of diagrams, and to 
sample (collections of) these diagrams by a Monte Carlo 
procedure. We split the Hamiltonian H of the impurity 
model into two parts, H\ and H 2 = H — Hi, and employ 
an interaction representation in which the time evolution 
of operators is given by Hi : 0{t) — e rHl Oe~ rHl . In this 
interaction representation, the partition function can be 
expressed as a time ordered exponential, which is then 
expanded into powers of H 2 , 



Z = Tr 



e -PBx Te -f£*rB*(T)] = g / dn--- dr n Tr\e-^-^(-H 2 )-.-e-^- T ^(-H 2 )e-^ H A. (12) 



Equation (12) represents the partition function as a sum over Monte Carlo configurations c = {n < . . . < r„}; n = 0, 
1, . . ., t, € [0, 0) with weight 



w r = Tr 



-(/3-T„)ff 



1 (-.Ha) • • ' e~ (T2 - ri)ffl {-H 2 )e- TlHl 



dT n . 



(13) 



Two types of expansions have been considered. In the 
weak-coupling approach |21j the partition function is ex- 
panded into powers of the interaction, H 2 — H; n t, while 
the time evolution between operators is given by the 
quadratic part of the Hamiltonian, Hi = H . The Monte 
Carlo configuration becomes a collection of interaction 
vertices on the imaginary time interval and the weight 
(131 evaluates to 



, weak 



= (-[/)" det 



G - \l_ 



dr r ' 



(14) 



Here (Go)y = Go(ri — Tj) is an n x n matrix whose el- 
ements are noninteracting Green functions evaluated at 
all time intervals defined by the vertex positions. Note 
that in the case of half filling of interest here, only even 
perturbation orders appear in the expansion. Away from 
half-filling, odd perturbation orders become relevant and 
Ising-type auxiliary fields must be introduced to avoid or 



reduce the sign problem. We will in this paper employ 
the continuous-time auxiliary field algorithm described 
in Ref. |24j . which for models with density-density inter- 
actions and an appropriate choice of parameters is equiv- 
alent to the weak-coupling algorithm [2"5] . 

The alternative approach is the hybridization expan- 
sion |22) where the partition function is expanded in pow- 
ers of the hybridization term, 



H 2 = y^(Vak*c 

aha 



f d 



+ h. C.), 



while the time evolution between operators is given by 
the impurity plus bath part of the Hamiltonian. This 
time evolution no longer couples the impurity and the 
bath. It therefore becomes possible to integrate out the 
bath degrees of freedom analytically to obtain 



= Z 



bath 



Tr 



x det M 1 {{n,ai}, {r„, a„}; {t{, a'J, . . . , {<, a' n }){drf 



(15) 



The configurations c are now collections of n time ar- 
guments Ti < ... < t„ corresponding to annihilation 
operators with flavor indices ai, . . . ,a„ and n time ar- 
guments t[ < ... < r' n corresponding to creation op- 
erators with flavor indices a[, . . . , a' n . The element i,j 



of the matrix M 1 is given by the hybridization function 
Pa'., a, { T i~ T j)i which is defined in terms of the hybridiza- 
tion parameters Vp' a and the bath energy levels [26] . 
In a model with density-density interactions only, one 
can separate the operators according to flavors, which 
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leads to the so-called segment representation [22 . This 
segment representation allows a simple and efficient eval- 
uation of the trace over the impurity states in Eq. (fl5 ). 



A. Implementation 

The implementation of the weak-coupling CT-QMC 
for the action ( 10 ) is straightforward. The non- 
interacting Green's function (111 is being Fourier- 
transformed and the resulting Go imn (r) inserted into 
Eq. fl4|). 



The implementation of the hybridization approach is 
more subtle, as - except in the equilibrium limit $ = 0, 
iip m = - the hybridization function F a >. a .(r[ — r,j) 



ing, because it is not directly related to the hopping am- 
plitudes V in the physical Hamiltonian 0. However, 
the hybridization function is implicitly defined by rewrit- 
ing the effective action (10 1 as [22] S e g — Sf + Si oc , 
with S F = -Eallo dr dr'd CT (r)F(r - t')4{t') and 



Si oc = -J Q dr (J2a £ ddld a -Udid[d t d l ). Consequently, 



the hybridization function can be constructed from (11) 
as 



F(-iu n ) = iu n -e d -G (iip m ,ioj n y 



(16) 
(17) 



which appears in the action (15 1 lacks a physical mean- After straightforward algebraic manipulation, we obtain 



FQw n ) = 



n [irsgn(w„ - § ip m )] 

a=±l 



[iw n - s d + a(i<p m - $)]^sgn(w„ - §<£ m ) - 



a=±l 



i<p m — $ 
2 



£ d + V s g n ( w n ~ f <Pm) 

a=±l 



(18) 



Note that the expression irsgn(w„ — aip m /2) emerges 
from imposing the wide-band limit for the leads. The 
hybridization approach is only able to cope with finite 
bands, because in the limit of infinitely wide bands of 
constant DOS, the expansion order diverges. The sgn- 
function must therefore be replaced by a sufficiently well- 
behaved function corresponding to a finite bandwidth 
and thus decaying rapidly enough for large frequencies 



by 



The high-frequency behavior of expression ( 18 ) is given 



F(h 



a=±l 



iT 

~2 



a 
2 1 



F(iu} n ) + - 



Cl 



Cl 



i(f m - $ 



(19) 



which means that the numerical evaluation of Eq. (17) 



requires some care. Conventionally, one regularizes the 
sum by analytically evaluating the dangerous parts and 
then numerically calculating the difference between the 
full function and the problematic parts, i.e. 



F(iw„) - F(iu n ) + t 



Cl 



10J, 



The leading order high-frequency tail Ci/(iw„) results 
in a constant shift — ci/2 in F(t), < t < /3. The first 



term 



a=±l 



iT 
~2 



<Pr. 



in the high-frequency expansion yields 

. _ T cos(^ m r/2) 



F(t) 



E 

71 — — OO 



F(iw n )e" 



j3 sin(7rr//3) 



(20) 



and diverges for r —¥ and r — > f3. These divergences 
are a direct consequence of the wide-band limit, i.e. we 
need to regularize them in order to be able to use the 
hybridization expansion algorithm. This regularization is 
introduced by cutting the divergences with a sufficiently 
large cutoff parameter F cut , i.e. we use 

F(t) = AF(t) - | + min (P(r), F cut ) . 

In practice, the value F cut = 10 4 was used. The contri- 
bution AF is Fourier transformed easily by accumulating 
the series numerically. 

Note that the term F has, besides the additional oscil- 



lations from the cosine modulation in Eq. ( 20 ) , the same 



structure as in the plain equilibrium Anderson model, 
where F(t) = ^(sin(T//3)) _1 . We will therefore illus- 
trate the properties of the quantity 



F(t) - F(t) - F(t) 
in the following section. 



(21) 
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FIG. 1: (color online) Imaginary-time data used as input for 
the CT-QMC solvers for different values of the imaginary volt- 
age ifi m - The upper panel shows the non-trivial contribution 
F(t), Eq. (21 1 to the hybridization function F(t), the lower 
panel shows the imaginary-time Green's function Go. Raising 
ifm leads to increasingly oscillating imaginary-time Green's 
functions and hybridization functions. The oscillations need 
to be resolved well by the QMC solver in order to guarantee 
an unbiased solution. As implied by Eq. ( 19 1 a strong neg- 
ative shift — Ci/2 occurs in the hybridization function when 
sweeping through the region tp m 3> The imaginary parts 
Im F(t) and lmGo(i~) are small and also show oscillations. 



B. Imaginary-Time Data 

Typical input data for both, the weak-coupling and 
the strong-coupling approach, are shown in Fig. [I] With 
increasing imaginary voltage tp m , oscillations with m 
nodes occur in both, the imaginary-time Green's func- 
tion and the hybridization function. Moreover, the shift 
( 19 1 grows quadratically, introducing a strong shift of the 



FIG. 2: (color online) Absolute values of the average sam- 
pling weight phases |(u> c /|w c |)| and |(wg/|wg|)| (Eqs. |14| 
and 



hybridization function towards negative values. 

The strongly oscillatory behavior for large tp m makes 
a correspondingly fine resolution of the imaginary-time 
interval necessary. In a standard Hirsch-Fye algorithm 
[27] . the interval [0, j3) has to be represented by a com- 
paratively small and fixed number of equidistant mesh 
points, i.e. these oscillations cannot be adequately re- 
solved. This limitation does not apply to CT-QMC, and 
it is hence the method of choice to access also large <p m . 



C. Phase Problem 

In contrast to the equilibrium case, complex sampling 
weights w c — e i7 |ui c | are obtained in both the weak- 
coupling and strong coupling formulation. As usual, one 
uses the modulus \w c \ of the weight to determine the 
acceptance probability, while the phase e n has to be 
treated as additional observable. Usually, such an ap- 



( |15[ )) of the weak-coupling (solid lines) and the strong- 
coupling (dashed lines) solver, respectively, as a function of 
the imaginary voltage. On a logarithmic scale, the average 
phase decays faster than linearly for the strong-coupling ap- 
proach when tp m is increased. No strong dependence on ip m 
is found for the weak-coupling algorithm. 



proach leads to a sign problem and severely limits the ap- 
plicability of the Monte-Carlo simulations. Therefore, we 
must anticipate a generalized sign problem, i.e. (e 17 ) — » 
exponentially or worse. The situation is especially prob- 
lematic for the hybridization expansion due to the ad- 
ditional shift (JT9J) towards negative values. Indeed, as 
illustrated in Fig. [5] the sign problem becomes increas- 
ingly severe with increasing imaginary voltage ip m , lim- 
iting this algorithm to small ip m . From Fig. [2] it also be- 
comes clear that the sign problem in the weak-coupling 
CT-QMC simulations is much milder and this approach 
allows us to simulate impurity models with large ip m . 

To demonstrate the quality of the imaginary-time data 
which can be obtained with the weak-coupling CT-QMC 
method, we show in Fig. ^ the imaginary part of 
the Matsubara axis self-energy computed for U = 10r, 
$ = 0.018r, T/T = 0.0098 and ip m = (m = 0), 
<p m /T = 1.23 (m = 10), (p m /T = 2.46 (m = 20), and 
Vm/r = 3.69 (m = 30). The equilibrium Kondo temper- 
ature for this parameter set is Tk/r sw 0.018 -C 1, i.e. we 
are reasonably deep in the Kondo regime of the Ander- 
son model. Moreover, the values for $ and T are such 
that T w Ik/ 2 and $ ss Xk, i.e. precisely in the parame- 
ter region which is hard or impossible to access for other 
methods. Even for large complex voltage the accuracy of 
the numerical data is very good (error bars on the order 
of the line width) for both small and large Matsubara 
frequencies. In contrast to the results presented in Ref. 
Il9j which are based on discrete-time Hirsch-Fye simula- 
tions, no discontinuities are observed for uj n w ±(p m /2 
in the CT-QMC data. We note, however, that a recent 
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FIG. 3: (color online) Imaginary part of the impurity self- 
energy obtained with the weak-coupling CT-QMC solver for 
V G = 0, U/F = 10, T/Y = 0.0098 and $/r = 0.018 . The 
equilibrium Kondo scale here is Tk/T ~ 0.018. We easily ob- 
tain high-quality data for all values m = 10, 20 and 30 of the 
complexified voltage, even in this most challenging parame- 
ter regime T K < I\ $ « T K and T « Ik/2. Each m-value 
was run on a single Intel Xeon(R) E5345 CPU for approx. 24 
hours, so the data were obtained with relatively moderate 
computational effort. 




FIG. 4: (color online) Geometric structure of the complex 
space carrying the two- variable Green's function G(z v ,z u ). 
Branch cuts occur for Ymz v = -Imz^, with j = ±1 (solid 
lines), for U — 0, but also at 7 = ±3 (dash-dotted lines), 
7 = ±5 (dashed lines), 7 = ±7 (dotted lines), and so on, for 
(7^0. Concentrating on the retarded sector of the Green's 
function, Imz w > 0, we introduce the cones bounded 
by the branch cuts with imaginary-part ratios and ^xr- 
Adding the real subspaces (Re z v , Re zj), the tubular cones 
are obtained as domains of holomorphy. 



preprint |28j reports a trick by use of which this issue 
could be resolved within the discrete-time formalism. 



IV. ANALYTIC CONTINUATION 



A. Analytic Structure 



As noted in Ref. [19] , at finite interaction, branch cuts 
occur for Imz w = ^Imz v (7 odd) in the complexified 
Green's function G(itp m — > z v ,iuj n — > z u ). Introducing 
the complex vector variable z= {z v ,z u ) we hence assume 
the Green's function to be holomorphic as a function of 
two complex variables in domains T C " :— M 2 +iC*, where 
for v G 2Z 



C' := 



sa > A 



<b< 



1 



are the cones emerging from the branch cut condition 
for positive (s = +1) or negative (s = —1) imaginary 
voltages (see illustration in Fig. [4]). Note that domains 
like T C " are well-known objects in the theory of functions 
of several complex variables and are called tubular cone 
domains. For a good introduction see, e. g., Ref. [5T] , 

In Ref. [H] this structure is described by the Cauchy 
representation 



S(icp m ,icj„) 



E 

7S2Z+1 



de- 



cr 7 (e) 



\UJ n 



(22) 



for the corresponding self-energy. However, Eq. (22) is 
only approximate, because the i</j m -dependence of the 
functions <7 7 (e) is not taken into account. Such a non- 
trivial dependence appears as a result of higher-order cor- 
rections in U. 

Let us start by discussing the analytically continued 
bare Green's function 



Co {ztp, z^) — 



a=±l 



r a /r 

f(z v -$) + ir(«)(z ,z w )' 



(23) 

with T 1 - ' (z v , z u ) := rsgn(Imz w — aImz v /2). The cor- 
responding geometric structure of the complex space is 
depicted in Fig. [4j the branch cuts given by the black 
lines 7 = ±1. Note that the Green's function does not 
vanish for all directions within a given T c » as |z| — > 00. 
On the other hand, ImGo(z) is at least bounded, and 
we assume that nonzero interactions do not alter this 
fundamental property. One can thus always find a con- 
stant c such that the imaginary part of the function 
f(z) := G(z) + ic is positive. Integral representations of 
the form J f(()K(z, £) dC = f(z) which are valid for the 
class of holomorphic functions with non-negative imagi- 
nary part also hold for G(z), since — ic • const (z) is also a 
function with non-negative imaginary part. This class of 
functions on tubular cone domains was extensively stud- 
ied by mathematicians. In Ref. [32], Vladimirov finds 
a generalization of Herglotz-Nevanlinna representations 
[3"3"] to such domains. See Appendix [A] for details. 

The validity of the imaginary-voltage formalism is 



presently based on the assumption of asymptotic conver- 
gence of the perturbation series in U. Thus, the influence 
of the branch cut between T C "+ 2 and T G « is expected to 
become negligible as v — > oo, i.e. all branch cuts with 
v > t'crit can be ignored. The maximal value v cr n may 
for example be estimated from the expansion order his- 
togram of the weak-coupling QMC simulation, since a 
given branch cut with index 7 = v + 1 is only established 
by diagrammatic contributions with order larger than a 
certain value n, which is roughly proportional to |7|. 
As stated in Ref. [TH] we are required to first take the 



Im z, 



limit 



$ and then 



i0 + . In our language, 



the spectral function is given by 



1 



A(uS) = lim lim lvaGt v \(z). 

7T v— >oo z— ►(<£>,«>) 



(24) 



Since branch cuts with index 7 > v CT n 
choose the domain T Ce with 



1 vanish we 



C £ := {(xi,x 2 ) G M 2 I x 2 > A -ex 2 < x x < ex 2 }, (25) 



and e 



for the analytic continuation of the inter- 



acting Green's function. This choice of domain is illus- 
trated in Fig. [5] In practice, the critical branch cut is yet 
chosen arbitrarily but to be small, see section IV B As 
shown in Appendix [X] the Poisson kernel representation 
resulting from Vladimirov's theorem is 



ImG{z)\ TCe = [ d 2 xV e (z-x) lim ImG(C) 

JR 2 



(26) 



with 



r. (27) 



7T 2 e ^iii ( X2 ~ M^iAO 2 + (2/2 - HVi/e)- 
where x and y are the real and imaginary parts of z. 

B. Maximum Entropy Method 

1. Single Analytic Continuation 

The numerical analytic continuation of imaginary-time 
quantum Monte Carlo data is a highly ill-posed problem. 
Even if the finite set of QMC data did not contain any 
stochastic noise there would exist an infinite-dimensional 
manifold of solutions to the integral equation associated 
with the continuation, i.e. the spectral representation 



G(iw„) 



de 



A(e) 



iu:,. 



K eq [A}(Lo n ) (28) 



for the conventional continuation problem. 

Hence, a regularization procedure picking a "most 
probable" solution is required. Typically, this is ap- 
proached with a Maximum Entropy Method (MEM), a 
rigorous framework rooted in Bayesian logic which can 
be understood as an automatic Ockham's Razor, in the 




FIG. 5: (color online) Sketch of the geometry of the two- 
dimensional analytic continuation problem. For the critical 
domain index i/ CI a the branch cut 7 cr i t + 2 = u CI i t + 1 (dotted 
line) is negligible, while the critical branch cuts ±7 cr it are not. 
The Green's function is therefore holomorphic in the cone 
domain T Ce bounded by the ratios ±Imz„ = — — Imz u =: 
dmz„, with C e given by Eq. ( |25[ ). Investigating the Green's 
function at the edge of this domain is compatible with the 
limiting procedure of taking z v — > $ and then z u — > w + 
i0 + for the spectral function A(lj). This is indicated by the 
bold dash-dotted arrow. Using the integral representation 
( |26[ ), a most likely limit of the Green's function at the edge, 
lim^_j. £ G(z), x 6 K 2 , will be inferred from the QMC data 
G(iip m , i<^n)| T c e in the domain using a Maximum Entropy 
Method (Sec. |IVB[ |. The spatial locations of the QMC data 
points in the domain are symbolized by the crosses. In the 
case of strong interaction, for small Matsubara frequencies we 
are limited to small values of (p m . 



sense of being "maximally noncommittal with regard 
to missing information" [34l436j . The spectral function 
A(uj) is interpreted as a probability distribution. A de- 
fault model D(lj) is introduced as a-priori information 
about the solution A{ui). Additional information, given 
by the measured imaginary-time data G(iw n ), is inferred 
through the kernel if e q[^4] in (28 1. If there is no addi- 
tional information the procedure will pick A(ui) — D(u>), 



in Bryan's MEM algorithm [37] . 
In practice, a functional 



Q[A] =x 2 [A\ -aS[A], a>0 



(29) 



is minimized in the space of candidate solutions for a 
given hyper-parameter a. The QMC data must be Gaus- 
sian distributed, such that the likelihood penalty ^[A] 
is given by 



1 N 



■K, 



eq[A]p)C pT f(G rj - 



■K eq [A] v ), (30) 



where G v are the measured mean real or imaginary parts 
of the imaginary-frequency Green's function G(iw n ), and 
Cp r j are the elements of the inverse covariance matrix. 
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The default model D(u>) is invoked through the entropy 



S[A] = J de 



A(e) - D(e) - A(e) lo{ 



Me) 
W) 



(31) 



For a detailed theoretical justihcation of this choice for 
the entropy see Ref. [36] . 

The easiest way of fixing the regularization parameter 
a is to employ the condition \ 2 ~ N (historic MEM). 
It is, however, more reasonable to calculate a posterior 
probability distribution Pr(A|a). Setting a to the max- 
imum of the posterior probability distribution is called 
classic MEM. Marginalizing a by choosing Pr(A|o;) as 
weights for A when integrating over a is empirically 
found to be most suitable and is also most justified from 
the theoretical point of view (Bryan's MEM). 



2. Double Analytic Continuation 

In order to adapt the above procedure to the double 
analytic continuation problem, a non-negative quantity 
has to be found which 

1. uniquely represents any possible function in the 
data range of interest - say T Ce - in order to define 
a x 2 f° r inference; 

2. easily allows calculating the non-equilibrium spec- 
tral function A(u>). 

We choose 



A(x) 



-- lim ImG«)| 

7T C— >x — 



(32) 



as such a representation, since due to the Kramers- 
Kronig relations and the validity of the representation 



(p6|), A yields a unique and simple representation of all 

The non-equilibrium spectral 
i($,w). 



possible functions G|yc, 
function is easily accessible, since A(uj) 
In the case of zero interaction, 



A (x) = --ha. V 

TV ' 



r Q /r 



a=±l 



x 2 -a{xi -$)/2-e d + ir 



(33) 



It is easy to verify that Aq (x) is a positive function with 
/ d 2 xA(x) = I if one constrains the xi-integration to an 
arbitrary finite interval of length I. This fact and the fact 
that A($,uj) = A(u) > do not imply A(xi,x 2 ) > in 
general. We however assume A(x\, x%) > and expect to 
obtain revealing signatures within the MEM, in case the 
real A is not positive definite for a given data set. Note 
that even in the presence of regions where A < 0, a MEM 
can be implemented, by identifying the nodes of A, as in 
the case of bosonic spectral functions. In general, posi- 
tivity may be enforced by adding a positive real constant 
b to the spectral function and adding a corresponding 



term to the image. As particular example for this proce- 
dure, we quote here the case of the Nambu off-diagonal 
Green's function Gi 2 , where the positivity is enforced as 
Gi 2 (r) + 6 / AuK{t,lo) = J c\ujK(t, u)(A a (u) + b) 02]. 

We hence choose ( p6| as a kernel function for the x 2 
functional and only take data in T e into account. The 
entropy expression (31) is adopted for a two-dimensional 
default model D(x). 



3. Implementation 

First note that since the input data for the Poisson 
kernel ( 26 1 are obtained from statistically independent 
QMC simulations, the covariance C in the \ 2 functional 
( 30 1 has a block-diagonal shape 



C 



c( m "-» +1 ) o 







V 



(34) 



C(mmax) J 



The submatrices G*™) are covariances for the subset of 
data G(i<p m , iw„) at a fixed tp m , estimated from the out- 
put of the corresponding equilibrium QMC simulation. 

Our implementation of the Maximum Entropy Method 
is based on Bryan's standard algorithm introduced in 
Ref. [57]. A singular value decomposition (SVD) of the 
kernel 



K : V- A -> V da 



K = VYAJ 1 



(35) 



is performed, with V, U T orthogonal, and the singular 
values 



£ = diag(cri,CT 2 , . . . ,<7 S ,0, . . . ,0), 



(36) 



ci > 02 > • • • > <y s > 0. Many important quantities may 
be reduced to the s-dimensional singular space Vs . Most 
notably, the (dim V^)-dimensional optimization problem 
given by 



Q[A] = min 



(37) 



may be solved within the singular space using Levenberg- 
Marquardt iterations. As s is comparably small after 
truncating the singular space with respect to the float- 
ing point precision of the singular values Cj (typically, 
s Rj 50), the algorithm is still sufficiently efficient, even 
though a two-dimensional frequency grid is required for 
the numerical resolution of A, and hence dim > 10 5 . 

The algorithm enables us to calculate several impor- 
tant data qualifiers and posterior probabilities and there- 
fore to classify both input data quality and candidate 
solutions. The posterior 



Pr(a|G) = Pr(a) / VA 



3 Q 



Z L Z s {a) ' 



(38) 
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Z s (a) 
-l 



/ VAe aS , and the 



with Z L = fV[KA]e~ x / 2 , 

Jeffreys prior [38] Pr(a) cx a -1 , is calculated using a 
Gaussian approximation for Q, centered around the so- 
lution A op t,a of Eq. (37). 

The usual procedures and strategies for data qualifica- 
tion and improvement of results as described in [35] are 
adopted: Assuming a flat prior Pr(D), the posterior for 
the default model 

Pr(£|G)cx J da J PiPr(a) ^ (39) 



had to be chosen in our case in order to take all relevant 
search directions in the A space into account. Quadruple 
precision floating point arithmetic was found to be un- 
necessary. For discretizing the A(x) function, logarithmic 
meshes for the x\ and X2 variables were used. Although 
A(x) does not decay for all directions as x — > oo, choos- 
ing a finite mesh and truncating the integrals was not 
found to be critical. 



4- Equilibrium 



is computed easily. Pi(D\G) serves as evidence for the 
quality of prior information when comparing within sets 
of default models for given QMC data. Whereas a pos- 
terior probability for the domain parameter e for given 
data and given default model, Pr(e\G, D), would be a 
sensible extension to the algorithm, we have not derived 
it yet. Useful ingredients might be found in the literature 
on blind deconvolution in signal processing, see |40) . In 
our implementation, a small enough e is chosen a priori. 

Picking appropriate data sets with well-estimated co- 
variance from the QMC output is also a non-trivial part 
of the problem. A good check is to determine the most 
probable mock error rescaling a where the covariance C 
is formally substituted by a 2 C. If the most probable a 
("merit"), i.e. the solution of 



-^classic 



+ N a = N 



(40) 



deviates from 1 by more than a few tens of percent, the 
input data are rejected Xciassic * s tne X 2 value of the 
classic MEM solution, the number of data points N, and 
the number of "good" data points N g = ^ - 
with Xi the eigenvalues of 



c+Ai 



A i3 = 



/J** 2 / 2 
dAidAj 



(41) 



In practice, a maximal Matsubara frequency n max com- 
patible with the error rescaling merit was determined, 
and all data ImG(i^ m , iuj n ) in T Cc , with n < n max 
were used for inference. Presumably, better data selec- 
tion strategies do exist. For example, using independent 
measurements for Re G and taking them into account by 
using a Schwarz representation (see Appendix |A| could 
yield better results. Furthermore, the largest Matsub- 
ara frequency index 7i max could be determined for each 
ip m individually. The latter appears to be necessary for 
non-equilibrium data. 

For the truncation of singular values, a threshold A was 
used, 



I Oi, if en > Acti max{Af, N}, 
^ 0, else (42) 



for an M by N kernel matrix. While for the conven- 
tional Wick rotation A ss 10~ 8 was sufficient, A « 10 -12 



As a test case we consider the equilibrium limit $ = 0. 
The data for (p m = can be analytically continued with 
the standard Wick rotation, using Eq. ( [28] ) and the stan- 
dard MEM. Figure [6] compares this ID spectral function 
to the result based on the 2D data set G(iip m ,iuj n ) and 
continued using the domain T B and the kernel function 
defined in Eq. (26). As default models for high temper- 
atures we use Lorentzians with variable width Tdefauit- 
They read 



£>(w) 



default 



i 



7TUj2 + r iefault 



for the ID continuation and 



D(x,u) = 



r. 



defau 



it fa) 



ttw 2 +r default (a;) 2 



(43) 



(44) 



for the 2D continuation, with Tdefauit^) = 
\/^dcfauit + x2 - ^ n annealing procedure in the tempera- 
ture was used for both, the ID and 2D data for invoking 
adequate prior information, i.e. we used the A solution of 
the next higher temperature as default model, starting 
with the Lorentzian at the highest temperature. This 
default model selection procedure appears not to have 
any strict Bayesian justification, however the physical 
argument is freezing out the high-frequency degrees of 
freedom and using present data for inferring low-energy 
details of the spectrum step by step [42]. A similar idea 
plays the key role in several modern renormalization 
group techniques. Note that Gaussian default models 
are not well-suited for our data, since the high-frequency 
tail in the wide-band limit is Lorentzian. This manifests 
itself quantitively in the following way: For the Gaussian 
default models we tested all had Pi(D\G) one order 
of magnitude lower than the Lorentzian ones. For 
both, the Gaussian and the Lorentzian, we can expect 
the quantity riefault /T to be > 1, due to the overall 
broadening introduced by a finite interaction U. 

Indeed, for the parameters U — 5T, Vq = 0, $ = 
shown in Fig. [6|jd), the (unnormalized) posterior prob- 
abilities Pr(D|G) and Pr(D\G) as a function of the pa- 
rameter Tdefauit are peaked at ss 2r for both, the ID and 
2D continuation procedures, respectively. These proba- 
bilities were calculated for j3T ~ 1QP. The most probable 
Tdefauit was chosen as default model. However, a strong 
dependence of the results on P e g was not observed. 
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FIG. 7: (color online) Lorentzian default model (441 with best 
Pr(Z)jG) for first annealing step in equilibrium. 



FIG. 6: Analytically continued data for equilibrium (<E> = 
0) obtained using the conventional Wick rotation (ID) and 
the unconventional two-variable continuation (2D), with U = 
5T, Vg = 0. The domain parameter for 2D continuation is e = 
j|. Subfigure (d) shows the posterior probabilities Pr(D\G) 
of the default models as a function of Tdefauit (Eqs. (431 and 

(PHI). 



A(Rez)T 
0.4 



The spectral functions shown in Fig. [6] were obtained 
for /3T = 5, 10, and 20. We chose e = ^ for the 2D 
domain, using n max = 10, 20, 40 for j3T — 5, 10, and 20, 
respectively. Note that due to the simple data selection 
strategy described in the previous section we only took 
into account data points with <p_ 2 < <p m < if2- Using 
a global n max , the estimate for the covariance subma- 
trix C'f" 1 ^ ) in Eq. (34 1 eventually becomes singular, even 
though C( m > with [ rn\ > 3 and uj n > w„ max could still be 
estimated for a limited set of Matsubara frequencies. We 
expect that using such additional, well-estimated 
might lead to more structured spectral functions. In 
practice, however, the merit a must yet be viewed as 
a rather crude measure of the quality of the covariance 
estimate. So for the purpose of both simplicity and re- 
producibility we used the stronger restriction. 

The Kondo temperature for U = 5T is Tk/T m 0.1, i.e. 
we can expect first signatures of strong coupling physics 
like Hubbard bands and a temperature dependent quasi- 
particle peak of reduced width in the spectra. Indeed 
both the ID and 2D MEM reproduces these features. 
More importantly, the overall shape of the spectra ob- 
tained agrees for all temperatures shown in Fig. [6ja-c). 
The results depend only slightly on the choice of Tdefauit 
for relevant values of Pr(D|G). Although the spectra 
inferred from the 2D procedure using our current imple- 
mentation appear to be less structured, the overall shape 
seems to be reconstructed quite well. For more serious 
calculations, the detailed high-frequency behavior (and 
behavior for large x) should be introduced with a more 
sophisticated default model, e.g. based on perturbation 
theory. 




Re ZhJ /T 



FIG. 8: (color online) MEM solution for A inferred from the 
QMC data at the lowest temperature, /3T — 20, for the equi- 
librium test case shown in Fig. [6] 



5. Inferred Representation 

Figure [7] shows the Lorentzian default model we used 
at the highest temperature in the annealing procedure, 
pT = 2. At the lowest temperature j3T = 20, the repre- 
sentation shown in Fig. [8] was obtained. The equilibrium 
spectral function shown in Fig. [6^c) is given by the cut 
A(& = 0,ui). Other values of Rez^, do not have any 
physical meaning. Note that certain structures appear 
in the inferred A(x, y) which vary as the domain param- 
eter e is changed: they occur for ~Rez v = zLKeez^. We 
interpret them as resulting from the properties of the 
kernel function discussed in Sec. IIV B 81 in combination 
with the MEM principle of only incorporating changes 
which are strongly supported by data. Also, at larger 
distance from the origin, discretization errors from the 
discretization of the double integral are most dominant 
for this most structured region of the kernel. At finite 
bias, the qualitative structure of the inferred representa- 
tion remains unchanged. 
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6. Finite Bias 

The rule of thumb n max ^- appeared to be a good 
choice for preparing the equilibrium QMC data for infer- 
ence. For $ > a first interesting observation is that at 
sufficiently low temperatures n max seems to be consider- 
ably smaller than 4^ . 

In fact, the simple data selection strategy yielding n max 
does not appear to produce a sufficiently informative data 
set to obtain quantitative agreement with for example 
RT-MC calculations [4 j. We observed this problem for 
f3T = 10 and the interaction strengths U = 4T and 
U = 6r and several values of the bias voltage $. On 
the other hand, by picking an n max for each ip rn sep- 
arately, we found larger sets of admissible input data, 
which tend to show a good agreement with RT-MC data 
for the current-voltage characteristics. While the pro- 
cedure is yet somewhat arbitrary, the following criteria 
were used to restrict the choices of data sets producing 
convergent MEM solutions: 

• ensure an error rescaling a ~ 1; 

• discard strongly oscillating solutions and solutions 
with obvious artifacts around uj ss 0; 

• discard solutions which strongly violate the physi- 
cal sum rule := J dui A(u) = 1. In many cases, 
too small values \\A\\ w 0.9 were obtained. Note 
that the MEM as we implemented it only has prior 
information about the value of the truncated dou- 
ble integral J J d 2 xA(x), because two-dimensional 
probability densities are considered when the en- 
tropy expression ( |31| is straightforwardly general- 
ized with respect to A; 

• use as many data points as possible, starting with 
small oj n , to maximize the amount of accessible in- 
formation. 

Note that the domain parameter e was, again, chosen 
somewhat arbitrarily: For U = 4T we only investigated 
I'max = 16, for U = 6r we picked v max = 20, with e = 
- — The dependence of the results on the particular 
choice of e was not studied systematically yet, but work 
along these lines is under way and the results will be 
presented elsewhere. The usual annealing procedure with 
temperatures (3T = 2,0T = 5, j3T = 10, where for f3T = 2 
the Lorentzian default models with Tdofauit = 1.5T (U — 
4T) and r de f au i t = 2.ir (U = 6r) were found to be most 
suitable based on the posterior Pr(D\G). 

The current J was computed using Meir and 
Wingreen's equation [25] 



J = 



Jmax / doj [f L {u) - f R (u)] A(u), (45) 



0.2 



0.15 
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_ _ U=4r, RT-QMC (voltage quench) 
■ - U=6r, RT-QMC (voltage quench) 
+ U=4r, Imaginary-Time QMC + 2D-MEM 
X U=6r, Imaginary-Time QMC + 2D-MEM 
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FIG. 9: (color online) Current-voltage characteristics ob- 
tained using the 2D MEM compared to RT-QMC g] data 
for indicated Coulomb interactions at temperature j3T — 10. 



Our experience up to now indicates that for too small 
sets of QMC data the method systematically underesti- 
mates the current, because Bryan's algorithm by conven- 
tion does not incorporate any changes to A D in case 
the data do not provide sufficient evidence for such mod- 
ifications. As a result, the current is too small, because 
in the vicinity ofw~0 the less structured default model 
obtained from the next higher temperature (initially the 
broad Lorentzian (44 )) is much flatter than the true solu- 
tion, which features a sharp Abrikosov-Suhl resonance in 
the relevant frequency range. Hence, the spectral func- 
tion obtained from the MEM has less spectral weight in 
the integration window in Eq. (451 than the true A(u>). 

Due to this trend of underestimation, in Fig. [9] we com- 
pare the largest values of the current compatible with the 
above-listed restrictions to data obtained using a recently 
developed RT-MC approach [3] . A generally good agree- 
ment is obtained. However, the data selection procedure 
is still too arbitrary to consider these results unbiased. 
Error bars are not available. If we only considered a 
fixed set of data G, the covariance Gov(A(x^), A{x^>)) 
would be estimated easily [37] . However, due to large off- 
diagonal terms, attempting to estimate an error bar for 
J is rather cumbersome. The $/r = 0.0625 run did not 
converge to a solution meeting our criteria for U = AY. 



with J, 



- £e 
max — ^ 



7. N on- Equilibrium Spectral Functions 

Spectra resulting from the procedure described above 
are shown in Fig. [10] These are the spectral functions 
used to compute the current in Fig. [9] While oscilla- 
tions appear, presumably due to the neglected error of 
the covariance estimate 39] , it is evident that the over- 
all spectral weight at small uj is larger for U = 4T than 
for U — 6r when <f> < 0.5r. This is consistent with the 
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FIG. 10: (color online) Spectra A(ui) = A($,u) used for the FIG. 11: (color online) Spectra A(oj) = A($,u) for U = 4F 
computation of the current shown in Fig. [9] as compared to fourth-order perturbation theory. 
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TABLE I: Norms of the spectral functions shown in Fig. |10| 



expectation that the quasi-particle resonance for U = 6T 
is already suppressed, because = 0.1 > Tr, whereas 

fi- 



ll 



we show a comparison 
= 4 and (3Y = 10 to 



T K for U = AY. In Fig _ 
of the spectral functions for U/Y 
the result obtained from fourth-order perturbation the- 
ory [3J. Based on the results presented in Ref. HI we 
expect that fourth-order perturbation theory is quite ac- 
curate at this interaction strength and temperature. Be- 
sides the unphysical oscillations in the MEM result and 
a bias towards the high-temperature default model, es- 
pecially for larger voltage biases, the agreement between 
the spectral functions, in particular the qualitative dis- 
tribution of the spectral weight, seems satisfactory. 

Table [i] shows the norm \\A\\ = J A(u>) dui for the func- 
tions presented in the figure. Obviously, the physical 
sum rule \\A\\ = 1 is not strictly obeyed, and there is 
a slight tendency towards too small norms whose origin 
is unclear but which appears to be consistent with the 
trend of current underestimation. Moreover, the selec- 
tion of data we chose at j3Y = 10 for U = 4Y and U = 6T 
is shown in tableland table |Hl[ respectively. The tables 
present the number N m « n max (m) — 2m je of Matsub- 
ara frequencies which are located within the cone domain 
T Ce for the chosen n max (m). We did not consider larger 
values of to, although at least m = ±4 yields further rel- 
evant information about A. For a test case the spectra 
did not show dramatic qualitative changes as additional 



TABLE II: Number N m of Matsubara frequencies taken into 
account for each value of m taken into account in the data 
selection at j3T — 10 for U/T = 4 and for the voltages plotted 
in Fig. H 



values at larger ip m were included, as long as the error 
scaling merit remained a ~ 1. However, the level of ar- 
bitrariness in the data selection would have been even 
larger, because of the corresponding additional n max pa- 
rameters. 

Obtaining reliable spectral functions at finite bias will 
obviously require more effort and we will briefly comment 
on possible avenues for this effort in the Conclusion. 



Kernel Structure 



We finish with some remarks about the structure of 



the kernel function ( 26 ) and its role in the continuation 



$/r 


N m =o 


N m =±i 


iVm=±2 


iVm=±3 


0.0625 


20 


11 
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1 


0.125 


21 


11 


6 


8 


0.25 


21 


11 


6 


3 


0.5 


20 


11 


6 


1 



TABLE III: Same as Table |T] but for U/T = 6. 
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problem. In the language of Bayesian inference the kernel 
function defines the information channel through which 
evidence about the shape of the representation function 
A(x) and thus also the physical spectral function A{uS) is 
extracted from the Monte Carlo data. 

For the information provided by a single data point, 
this channel results in vague (strong) evidence for 
changes in a given compact region R C V^, see Eq. (35 1, 



depending on whether the subset of column vectors Ui 
of U spanning R is associated with small (large) singu- 
lar values <7i, and a small (large) overlap of the column 
vectors Vi of V with the data point. For this reason, 
very small singular values yield irrelevant components of 
the channel and are therefore projected out in Bryan's 
algorithm by introducing the threshold A, Eq. (42). 

We can neither perform the SVD analytically, nor 
can we analytically take into account structural changes 
which occur when rotating the basis of Vdata to the eigen- 
basis of the covariance matrix C in order to consider sta- 
tistically independent data. We can however consider 
values of the kernel in Vfi for a given data point, as- 
suming it to be uncorrelated with other data points so 
that it may be investigated separately. Within our QMC 
implementation, experience shows that correlations be- 
tween Matsubara frequencies ui n , uj n > are monotonically 
decreasing as a function of distance \uj n i — u) n \, though 
very slowly. 

Let us first consider a single uncorrelated imaginary 
part of a Green's function at Matsubara frequency u> n 
in the standard Wick rotation problem. The spectral 
function A(e) is inferred through the Lorentzian-shaped 
kernel ((281), 



lmK eq [A(e)](u n 



(46) 



For all uj n the kernel (28) is centered around e = and 
higher frequencies are associated to larger values of the 
kernel as the width given by w n is increased. As com- 
pared to e ps the values of the kernel at large frequencies 
are still small. We can therefore expect large singular val- 
ues and thus relevant components of the kernel to be as- 
sociated with small frequencies only. This is in agreement 
with the well-known observation that high-frequency in- 
formation about the spectral function is better put into 
the default model as prior knowledge and a good resolu- 
tion is obtained for the - fortunately most interesting - 
low-frequency region. 

In the case of our two-dimensional continuation the 
situation is quite similar. For given data Im G(i(p m , iuj n ) 
the Poisson kernel in Eq. ( 26 ) is 



n 



LO n - flipm/s 



t^ x (x 2 - L-txi/e) 2 + (uj n - ^ip m /e) 2 



It is the product of two Lorentzians. In analogy to the 
argument given above one may expect the best resolution 
for data i(x< best )) with 
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FIG. 12: (color online) Expansion order histogram obtained 
using the weak-coupling solver for K = —/3U/4 + 1, which 
suppresses the odd perturbation orders 4]. The results are 
for U = IT, PT = 51.2, and <^ m = 2.46r for indicated voltages 
$. The average order decreases as $ is increased. A similar 
behavior is obtained for different values of (p m and U. 



This does not depend on the physical voltage $, except 
that the critical branch cut index 7 cr j t appears to be de- 
creasing as a function of This can be estimated from 
the expansion order histogram for the example shown in 
Fig. [12] Consequently, the domain parameter e could 
presumably be raised as <& be increased. However, in the 
limit of very large voltages, especially the low-frequency 
region of the physical spectrum A(uj) — A(§,uS) is not 
expected to be in the best resolvable region ( 47 ) . 



(best) 



■ (best) / n (best) (best) 

±x^ /£ and x\ ,x\ 



0. 



(47) 



Thus, the approach based on a representation of data 
in T Cs appears to be limited to relatively small volt- 
ages. Note that, since $ w Tk is the most interesting 
parameter regime, this is presumably no serious draw- 
back. However, identifying subtle details in the range 
— 4>/e <C u) <C $/e may require more care than the case 
uj fa for the standard Wick rotation. Fixing the X2 and 
x\ variables in the kernel and analyzing the dependence 
as a function of the data coordinates (f m , ui n we similarly 
find that large values of the kernel are found in the vicin- 
ity of the domain boundary, i.e. for (m, n) pairs close to 
the cone boundary, ui n ps ±(p m /s, with ip m , u> n not being 
too large. Hence, data close to the boundary provide the 
most relevant information. This appears to explain the 
importance of an m-dependent n max in our computation 
of non-equilibrium spectra. 



V. CONCLUSION AND PERSPECTIVE 

The imaginary time formulation for steady state trans- 
port in strongly correlated quantum impurity systems 
proposed by Han and Heary is based on the solution of 
a family of quantum impurity models subject to com- 
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plex voltages, and a subsequent double analytical con- 
tinuation with respect to frequency and voltage. A main 
purpose of the study presented in this paper was to in- 
vestigate to what extent an unbiased, numerical imple- 
mentation of this approach is feasible and whether or not 
it yields physically plausible results. 

To solve the impurity problem we employed two re- 
cently developed continuous-time impurity solvers. The 
hybridization expansion approach was found to be un- 
suitable in the case of large complex voltages, due to 
a serious sign problem resulting from the shift of the 
hybridization function to negative values. The weak- 
coupling approach, on the other hand, works well for 
small and large ip m . Even though the non-interacting 
Green's function Go becomes complex and oscillating, 
the resulting sign problem is mild, enabling us to obtain 
highly accurate, unbiased imaginary-frequency data for 
all relevant complex voltages. This part of the problem 
can be considered as solved, leaving us with the double 
analytical continuation problem. 

A main result of this work is the derivation of an an- 
alytical expression of the kernel (Eq. ( [27] ) ) for the ana- 
lytical continuation procedure. This kernel is consistent 
with the analytical structure (branch cuts) of the theory 
and maps a function of two variables, A{x\ 1 X2) 1 to the 
interacting Green's function in a tubular cone domain 
of the complex voltage and frequency space. The phys- 
ical spectral function for a dot under voltage bias $ is 
obtained as A{oo) = A(Q,uj). 

We have implemented and tested an analytical continu- 
ation procedure based on the Maximum Entropy Method 
and our proposed kernel. We want to emphasize that 
both the data selection procedure and the estimate of 
the covariance entering into the maximum entropy em- 
ployed for $ > are at this point still rather rudimen- 
tary and leave room for improvement. Our results for 
the non-equilibrium case should therefore be viewed as 
preliminary and illustrate the presently most plausible 
spectral functions and currents which can be obtained 
using our current implementation. 

Nevertheless, taking into account the obvious chal- 
lenges inherent in a double analytical continuation proce- 
dure, we find physically reasonable spectral functions for 
the interacting equilibrium model and, to a lesser extent, 
also under finite bias. A comparison of the spectral func- 
tions with fourth-order perturbation theory shows that 
the approach is able to reproduce the correct trends, al- 
beit the strong oscillations resulting from the maximum 
entropy approach render a detailed comparison mean- 
ingless. On the other hand, the current calculated using 
these spectral functions is in fair agreement with recent 
results from a real-time Monte-Carlo approach. 

We hope that further improvements in data selection 
strategies, a better understanding of the precise behav- 
ior of the Green's function across the branch cuts, im- 
proved default model functions and, very importantly, 
the inclusion of the sum rules into the maximum entropy 
algorithm will eventually enable us to obtain more ac- 



curate results and turn the combination of Monte-Carlo 
and double analytical continuation into a reliable tool for 
the study of steady-state properties of quantum impurity 
systems using Han and Heary's formalism. 
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Appendix A: Derivation of the Kernel 



Based on the argument given in section IV A we re- 
strict ourselves to the class of functions with positive 
imaginary part in the domain T e , typically denoted 
as H + {T Ce ) in the mathematical literature. For a good 
overview of the concepts and terminologies used in the 
mathematical context see Ref. [30] and the first volume of 
Ref. [31] , Vladimirov found the following generalization 
of Herglotz-Nevanlinna representations to several com- 
plex variables [33] . It is essentially [3T] the 

Theorem. (Vladimirov, 1978/79) The following con- 
ditions for a function / e H + (T C ) are equivalent for a 
cone C C R m and fi(x) := Im/(x): 

1. The Poisson integral Pc[d/i] is pluriharmonic in 
T c. 

2. the function Im/(z), z = x + iy G T c , is repre- 
sented by the Poisson formula 

lmf(z) = P c [dfi](z) + (a,y), (Al) 

for some a € C* , where C* is the dual cone of C; 

3. for all z Q 6 T c , under the assumption that C is 
regular, the Schwarz representation 



f(z)=i S c (z-t_,z° -t)dfi(t) 
+ (a, z)+b 
holds, with b = b{z°) = Re/(z°) - (a,x°). 



(A2) 



□ 



Let us introduce the relevant mathematical terminology. 
A cone C C R m with vertex at zero is defined [3D] by 
the property that y e C =>• VA > : \y E C. Its 
dual cone C* := {£ e R m | Va; £ C : > 0}. Here, 

Pc[d/j,](z) = J Rm d m x \i{x)Vc{z — x) with the Poisson 
kernel 



V c {z) 



(2^)™/C c (2hy)' 



z = x + iy 



(A3) 
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and the Cauchy kernel 



Obviously, 



ICc(z) 



d m £e 



(A4) 



c* 



We will not explicitly use the Schwarz kernel S, the 
reader may find it in Ref. 31 . A holomorphic map- 
ping is said to be biholomorphic iff it is one-to-one. Two 
domains G, G are biholomorphically equivalent iff a bi- 
holomorphic mapping G — > G exists. For the concept of 
pluriharmonicity see introductory volumes of Ref. |31j . 
In the case of T Ce we rewrite Eq. ( 25 ) as 



C e = (J {(xi,x 2 ) S R 2 |x 2 > OAxi = Ax 2 }. (A5) 

Ae(-e,e) 

Hence, the dual cone 

C* = f) e R 2 |Vx 2 > : £ x Ax 2 +&x 2 > 0} 

= {(Ci,6)€K 2 |f 2 >0A6€ 



Evaluating the integrals / d m £ = J°° d£ 2 /_ 



in (A4) yields 



€2/e 



da 



KcM)=-- n 



(A6) 



Eq. (27) follows immediately from the definition rtA3 




In order to prove the validity of the representation ( |26[ ) 
based on Vladimirov's theorem, we first determine a = 
due to the boundedness of the Green's function. Now 
we need to show that the Poisson integral Pc c [ d/i ] with 
respect to the measure /i(x) — Im/(x) is pluriharmonic 
for all functions / e H + (T Cc ). Note that for the Tri- 
dimensional octant 



= {( 



Xi, 



> 0} (A7) 



it was proven [31] |4T] that the Poisson kernel V n (m.) is 

pluriharmonic for all functions / € H + {T C + ). For- 
tunately, as we restrict ourselves to m = 2 in our ap- 
plication, all tubular cone domains are known to be bi- 
holomorphically equivalent - they are simply connected 
through linear transformations. 

To see the advantage more explicitly, we introduce the 



biholomorphism M : T 
operation 

M(I) := M ■ I ' 



T c ' given by the linear 



(l + £ 2 )/2 



\-l (e- 



L )/2 y 



(A8) 



M~ 



(e + e^)/2 -(l+ £ 2 )/2 x 



l 



(A9) 



We explicitly show that the kernel representation ( 26 ) for 
a function f(z) € H + (T Ce ) may also be derived by apply- 
ing the corresponding Poisson kernel V r (i) for the tubular 



octant to the corresponding function f(z) := f(Mz) £ 
H+(T + ) and transforming back to T + . Since the 
representation for / is valid, we will have shown explic- 
itly that ([26]) is valid for all / e H + (T C -). 
For this purpose it suffices to show that 



JC Ce (z) =/C„ (2) (M- 1 z), 



(A10) 



because then Vc e (z — x) = ~P C (V (M 1 z — M 1 x) and 

therefore V c U - Mi) = V n ^)(M~ l z - x). We in- 

troduced the integration variables x and x of the Pois- 
son integrals Pnldfi], (x(x) = Im/(x) and P„(2>[d/2], 

+ 

/u(x) = /(x), respectively. Since det M = 1, transform- 



ing x — >• x in P„{2) then yields (26). 

°+ 1 — ■ 

With a similar procedure as for JCc e it is straightfor- 
ward to show that 



/C„(2) (2) 



1 



5 e t c + . 



(All) 



To finish the argument we verify that Eq. ( A10) holds by 
inserting 



/C r<2 )(M- 1 z) =(l + £ 2 ) 



e + e 



1 1+e 2 
— Z\ 



-Z'2 



■ (zi + ez 2 ) 1 . 

Representations for any tubular cone domains in C 2 are 
similarly related due to the biholomorphic equivalence. 
In particular, valid representations for T C " are obtained 
easily. For example, the Poisson kernel with respect to 
T c i reads 



n*) = i n 



2/2 - + m)2/i/ 2 



^ ui±x to - ^i) 2 + (y 2 - ^yi) 2 (A12) 



and could in principle be used for an enhanced contin- 
uation procedure invoking data from all sectors of the 
complex space. 
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